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We propose a simple method that allows, in one dimension, to solve exactly a wide class of clas- 
sical stochastic many-body systems far from equilibrium. For the sake of illustration and without 
loss of generality, we focus on a model that describes the asymmetric diffusion of hard core parti- 
cles in the presence of an external source and instantaneous annihilation. Starting from a Master 
equation formulation of the problem we show that the density and multi-point correlation functions 
obey a closed set of integro-differential equations which in turn can be solved numerically and/or 
analytically. 



Recently, nonequilibrium statistical mechanics in low 
dimensional systems has attracted much interest [jp. 
Whereas classical statistical mechanics has a well estab- 
lished body of concepts and tools to deal with systems 
at or near thermodynamic equilibrium, the situation is 
less satisfactory for far from equilibrium systems. The 
theoretical understanding of classical stochastic many- 
body systems where both purely diffusive motion (due 
for example to static disorder, phonons, etc.) and reac- 
tions (for example fusion, annihilation, coagulation, etc., 
where the number of particles is not conserved) take place 
simultaneously is relevant to a wide class of phenomena 
in physics, biology, economics, a.s.o. 
In view of this, we have chosen to investigate an im- 
portant class of models describing diffusion-limited re- 
actions (for recent reviews see 0). While it has long 
been recognized that the mean field approach is not ap- 
plicable in one and two dimensions, it is only recently 
that much progress has been achieved in the description 
of the anomalous kinetics of low dimensional systems. 
Some open problems in this field are reminiscent of those 
studied intensively in the field of strongly correlated low- 
dimensional systems. Hence, one might hope to use with 
advantage the powerful tools devised in the latter. 
The last decade has witnessed new developments in the 
Master equation approach. The Markov generator of the 
stochastic time evolution for the probability distribution 
can be interpreted as the Hamiltonian of quantum spin 
chains (with quantum group symmetries) or alternatively 
of fermion or boson many-body systems with interaction 
terms. The dynamics is then coded in an imaginary time 
Schrodinger equation or alternatively formulated as an 
Euclidean field theory. The most formal approach ex- 
presses the Hamiltonian as a sum of generators of an 
Hecke Algebra and, through Baxterization, associates an 
integrable vertex model The free fermion models 

considered by Grynberg et al. § , Schiitz and others (see 
[pj and references therein) can be viewed as special cases 
of the above with a grading of the Yang-Baxter algebra. 
As is known in the field of integrable models, it is in 
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general extremely difficult, if not impossible, to extract 
the correlations from the Bethe "wave functions", even 
in their algebraic form (see for example Q). This is un- 
fortunate because conceptually as well as experimentally 
(see Q and references therein), the correlations coded 
in these "wave functions" are the most relevant quan- 
tities. We should mention at this point however that 
there is one successful application of ideas developed in 
the field of strongly correlated systems, namely the ex- 
tension of the numerical Density Matrix Renormalization 
Group to nonequilibrium statistical mechanics in one di- 
mension [Q. 

On the analytical front, a powerful method is that devel- 
oped by Cardy and collaborators ||. Assuming that the 
system is in a low-density regime, these authors rewrite 
the Markov generator of the stochastic time evolution as 
a bosonic Hamiltonian which in turn can be treated by 
field theory and renormalization group techniques in ar- 
bitrary dimension. For certain models, it is even possible 
to evaluate exactly the beta functions for the coupling 
constants H . This is a first indication that the underly- 
ing mathematical structure is substantially simpler than 
that of the quantum many-body problem. To our knowl- 
edge, this observation as well as its implications do not 
seem to have drawn much attention in the literature. 
Despite the progress achieved in the last decade, there 
are still many open or poorly understood problems, even 
in one dimension: the multi-species and the epidemic 
models, the propagation of catalytic shock waves, etc. In 
response to these challenges, we were lead to consider a 
fcrmionic field theory formulation of the Master equation 
for dense systems || . While inquiring about the possible 
application of the renormalization group to the fermion 
field theory, we realized that it is possible to express the 
density and correlation functions as solutions of a closed 
set of integro-differential equations. These in turn can be 
solved exactly numerically and/or analytically. This is an 
unexpected and remarkable result, especially in view of 
the large amount of work that has been performed in the 
field using sophisticated methods. Our approach in con- 
trast is simple and elementary. In this letter we would 
like to report on this rather unanticipated observation. 
For illustrative purposes, we focus here on an extended 
Lushnikov model (l(J that describes the (asymmetric) 
diffusion of particles in the presence of bilinear source 
and annihilation. Certainly, such a model does not con- 
stitute a protoptype for classical stochastic many-body 
systems as does the Ising model for phase transitions: 
the variety of phenomena to be described in nonequilib- 
rium statistical mechanics is much too rich to be coded in 
a single model. Though we are perfectly aware of the nu- 
merous potential applications that our method provides, 
we will refrain from mentioning them here due to a lack 
of space and refer the reader to elsewhere (e.g. 
Consider a lattice of L (L/2 even) sites with periodic 
boundary conditions on which classical particles with 
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hard core can diffuse with rates h to the left and h! 
to the right, respectively. On encounter, pairs of par- 
ticles annihilate with rate e'. A source injects parti- 
cle pairs with rate e. Rewriting the Master equation 
as a Schrodinger equation in imaginary time, we have 



d\P(t))/dt = -H\P(t)), with (a 



m-\-L 



i) 



H = -J2 { ha m+i a m + h'a+a m+1 + e'a m a m+1 + e(a+a+ +1 - 1) 

m=l L 

+ (2e-h- h')<r+a m - (e' + e - h - h')cr+a m a+ +1 a m+1 \ (1) 



where er^ denote the standard Pauli matrices. The Dirac 
ket \P(t)} = J2 V •f»?(*)l r /)> w here r\ denotes a random con- 
figuration in Fock space, encodes the probabilities of the 
various configurations (see ||). The above dynamics is 
supplemented with the initial conditions |P(0)) = \p) 
(see below). It is convenient to introduce the left vac- 
uum (x\ — Ylnivl- Conservation of probability implies 
then (x\H — 0, while expectation values are linear in the 
kets (A(t)) = (x\A\P(t)). The models that have been 
solved previously correspond to e' + e = h + b! (see e.g. 
[||,||J|Jl^]), i.e., to quadratic hamiltonians. Below, we 
show that the Hamiltonian H can be solved exactly for 
7 = e? + e — (h + h') ^ by elementary means, i.e., that 
the density and two-point (multi-point) correlation func- 
tion^) can be computed exactly even though the hamil- 
tonian is non-quadratic. 

Performing a Jordan- Wigner transformation, and subse- 
quently a Fourier tranformation, leads to 

H = [ UJ (l) a l a q + ti(q) a - q a -q + 2sinq (e'a q a- q + ea^ q a q 



q>0 



eL 



1^2(1 - cosq)a f q a q + - ^ cos{qi ~ q2)a\ l a q2 a\ 3 a qi+q3 - q2 , (2) 

9 91,92,93 



where w(q) = a — bcosq + i(h' — h)smq, with a = 
e' — e, b — e' + e. Following reference 0, we consider 
a translationally invariant and uncorrelated initial state 
with an even number N of particles of density p — N/L, 



\P) 



2(1-^ 
1 + (1 - 2 P y 



m 

q>0 



P 2 cotq)al q al)\0) , 



(3) 



where p = p/(l — p). 

Using the 
left vacuum ((x\ a l = 



properties of the 

., ; — cot(|)(x|a_q, see [||), and with 
the notation (ci- q a q i)(t) = (x\a- q a q > e~ \p) , the den- 
sity and the two-point density correlation function can 
be written as 

p(t) = (ala m )(t) = ^Y. ei{q ~ q ' )mc <\)^\ a -^' e ' m \p) ( 4 ) 
9,9' 
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»>(*) = p E 



im(q— q' )-\-in(k— k' ) 



k,k ; ,q,q ; 



{a^ q a q >)(t)S k ',k ~ {a q >a k i){t)8 q - k 



+ cot(-)(a_ fc a_ g aya fe /)(i) 



cot(|) (5) 



Futhermore, it can be shown that 111 111 



[i 2 cot(f) 



(x|a- 9 V e I/ 5 ) =3(?>i){fl-?V>(°) =g ^ g, ^ l + /t 2 cot( a ) 



2^ X , 



(6) 



where g(g, t) is a function (even in q) that is a solution of a 
non-linear integro-differential equation (see below), with 
g(q, t = 0) = 1, i < 0) = and g(q, t) = g(q + 2tt, t). 
This remarkable property can be generalized to higher 
order correlation functions. Due to the conservation of 
probability, an analog of the Wick theorem applies and 
all multipoint correlation functions can be computed. For 
simplicity's sake we consider in the following a lattice 
initially full, i.e., p = 1 (for the general case see 

0) and 

obtain the integro-differential equation for g(q, t) 
d 

q19(q, t) = 2e(l + cos q) - 2(6 - j)g{q, t)+2(a--y) cos q g(q, t) 



/7T 
-TT 



27 

dk(l — cos k)g(k, t) g(q,t) I dk cos k g(k, t) 



27 r 

H g(q,t)cosq / dkg(k,t) 

with the initial condition g(q, 0) = 1. This equation 
with the appropriate boundary conditions fully solves the 
problem of the time evolution of the density and corre- 
lation function. Indeed, introducing the help functions 
F n (t) — (l/27r) f_ dq g(q,t) cos(nq) (for n integer ), the 
density can be written as 



(7) 



p(t) = F (t), 



(8) 



while it is convenient to introduce the function Q (r, t) = 
1 — 2 J2i< n <r ^2n-r-i(t) to express the (connected) two- 
point correlation function as 



C r {t) 



1 r 
4 



(9) 



g(r- l,t)G{r + l,t)-G 2 {r, t) 
The help functions F n (t) have the integral representation 
F n (t) = e-«^I n (K (t)) 

+ f dr(e + 7 {F x (i - r) - F (t - r)}) e ^(*-)-^(t) j 2 I n (R (t) - K (t - r)) 



in+l 



(K {t)-K {t 



_(K (t) - K Q (t - t) 



(10) 



where I n (z) denote the modified Bessel functions of 
imaginary argument, and Ko(t) = 2(a — 7)^ + 
4jJ*dsF (s), K^t) = 2(6-7)t + 47/ *daFi(«). It 
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is worth noting that Fo(t) and Fi(t) are obtained self- 
consistently while the higher order functions F n (t) (n > 
2) follow from simple integration. Notice that at 7 = 0, 
we recover the results obtained in §J|]|,|,|lO|] . We have 
solved the integral equations (10) and obtained the F n (t) 
numerically as well as analytically in the asymptotic 
regime. In the long time regime, the asymptotes follow 
as 

i) In the absence of source (ht, h't, e'f > 1 , £ = 0) 

1 7 

P{t) ~ y/4n(h + h')t ~ vre'(/i + h')t (U) 
-1 r 87 "1 1 

Cr{t) ~ 4n(h + h')t + V r + VI (4n(h + h>)t)V* {U) 

The result for the density agrees with that of reference 
[Oj. The numerical computation of p(t) is shown in 
Figfi 

ii) In the presence of source (ht,h't,\a\t,bt ^S> 1), two 
regimes have to be distinguished, namely 

1) If K (t) > 0, then, 

/x „ / s exp(— 4(e — 7p(oo) 2 )t) „. 
p(t) - p(oo) ~ -C r {t) ~ i y y t3/ P ' ' ' (13) 

A numerical illustration of this case is shown in Fig.|^ 

2) If K [t) < 0, then, 



exp 



(- 2[(b - 7) + 2 7 Fi(oo) - \a- 7+ 2 7 p(oo)|]i) 



p(t) - p((X>) ~ Cr(t) * '- (14) 

where p(oo) = 1/(1 + yj e'/ e ) i s the stationary density 
which coincides with the (free fermion) 7 = case, and 

Fx (00) = V«7(V^+^) 2 - 

The results jll|)-(14) indicate that the 7 term is irrele- 



vant in the long time limit in the sense that it does not 
renormalize the exponent of the algebraic (in the absence 
of source) or of the subdominant-algebraic (in the pres- 
ence of source) decay. Such a behavior had been observed 
numerically (in the absence of source) in reference [ fl2"| . 
The correlation times are however renormalized. 



The extension of our method to deal with epidemic 
models and many species is indeed possible. Local 
sources or sinks as well as boundary (initial state) terms 
associated with local operators in the hamiltonian formu- 
lation can also be dealt with whenever the Jordan- Wigner 
phase can be made to cancel. A physical illustration of 
this is realized in a model of photoluminescent diodes 
with injection currents at the boundaries. These and 
other physically interesting applications of the method 
proposed in this paper will be discussed elsewhere 111]] . 
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FIG. 1. Density p(t) versus time for e = 0, e' = 1.5 
and 7 = 0.5; (upper curve) 7 = (middle curve: 
free-fermions); 7 = —0.5 (lower curve). 
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FIG. 2. The logarithm of the density i i/2 \p(t) - p(oo)] 
(stars, p(oo) = 1/3) and of the two-point correlation function 
-t 3/2 C r (t), 1=1 (circles), r=2 (triangles) for e = 0.5 ,e = 2, 
and 7 = 1 versus time. We have included the shifted asymp- 
tote const — 14t/9 as a guide to the eye (thin line). 
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